Bayesian brain mapping BBM is a technique for producing individualized functional brain topographic maps from existing group-level network maps. Group-level network maps are commonly group ICA maps or group average parcellations, but other types of networks maps such as non-negative matrix factorization and PROFUMO can be used. In the case of ICA, BBM is known as template ICA (CITE Mejia et al 2020). BBM is a hierarchical source separation model, with priors on the spatial topography and, optionally, on the functional connectivity. The priors are estimated a-priori, so the model can be fit to individual-level fMRI data. The noise-reduction properties of the population-derived priors result highly accurate and reliable individual network topography maps and the functional connectivity between them. Importantly, the subject-level network maps are matched to the corresponding group networks from the template (i.e. parcellations or group ICA maps). Because BBM is applicable to individual-level analysis, it is computationally convenient and has potential clinical utility.
Once a set of group-level network maps has been chosen, there are two
steps to performing BBM. Both are implemented in the
BayesBrainMap R package.
Step 1. Training a population-derived prior (mean and variance) for each network using a training dataset.
Step 2. Fitting the Bayesian Brain Mapping model to fMRI data from an individual subject, using those priors.
Here, we perform Step 1 using data from the Human Connectome Project (HCP). Specifically, we train CIFTI-format Bayesian brain mapping priors using a variety of ICA and parcellation-based templates, listed below. The population-derived priors described here are available for use for individual-level Bayesian brain mapping. For analysis of individuals from other populations, it is often desirable to train the prior on a set of training subjects representative of that population. To facilitate this, we also provide and describe the code used to produce the HCP-derived priors, so that this workflow can be easily reproduced in other datasets. Finally, we illustrate the use of
For the choice of group-level network maps, we provide several options:
Group ICA maps from the HCP at resolutions from 15 to 50 (CITE)
The 17 Yeo networks (CITE)
To reproduce this workflow, first clone the repository to your local machine or cluster:
# git clone https://github.com/mandymejia/BayesianBrainMapping-Templates.git
# cd BayesianBrainMapping-TemplatesNext, download the required data/ folder from the
following OSF link:
https://osf.io/n3wk5/?view_only=0d95b31090a245eb9ef51fe262be60ef
Once downloaded, unzip the folder and place the data/
directory inside your cloned GitHub repository, so the folder structure
looks like this:
BayesianBrainMapping-Templates/ ├── data/ │ ├── template_rds/ │ └── parcellations_plots/ │ └── ... ├── src/ │ ├── 0_setup.R │ └── 1_fd_time_filtering.R | └── ... ├── BayesianBrainMapping-Templates.Rmd └── ...
This section initializes the environment by loading required packages, setting analysis parameters, and defining directory paths.
Important: Before running the workflow, you must
review 0_setup.R and install any necessary packages, ensure
you have an installation of Connectome Workbench, and update the
following variables to match your local or cluster environment:
dir_project (path to the main project folder that
contains all code and data (i.e., your local clone of
BayesianBrainMapping-Templates))
dir_data (path to the data/ folder
inside the project directory, downloaded from OSF in the previous
step)
dir_results (path to the folder where all output
files are saved)
dir_personal (path to your personal user directory
on the cluster)
HCP_restricted_fname (path to the restricted HCP CSV
if you have access to it)
wb_path (location of the CIFTI Workbench on your
system)
github_repo_dir <- getwd()
src_dir <- file.path(github_repo_dir, "src")
source(file.path(src_dir, "0_setup.R"))## Using this Workbench path: '/Users/nohelia/Downloads/workbench/bin_macosxub/wb_command'.Before estimating the BBM priors, we first select a high-quality, balanced subject sample to ensure reliable, representative priors. Starting from the full HCP sample of 1206, we apply the following filtering steps:
We begin by filtering subjects based on the fMRI scan duration after
motion scrubbing For each subject, and for each session
(REST1, REST2) and encoding direction
(LR, RL), we compute framewise displacement
(FD) using the fMRIscrub package. We use a lagged and
filtered version of FD (CITE Pham Less is More and Power/Fair refs
therein) appropriate for multiband data. FD is calculated from the
Movement_Regressors.txt file available in the HCP data for
each subject, encoding and session.
A volume is considered valid if it passes an FD threshold, and a subject is retained only if both sessions in both encodings have at least 10 minutes (600 seconds) of valid data.
The final subject list includes only those who passed the filtering
criteria in both LR and RL encodings and in both visits (REST1 and
REST2). This list is referred to as the combined list and
is the one used throughout this project.
# This script filters subjects based on motion using framewise displacement (FD) from fMRIscrub.
# For each subject, encoding (LR/RL), and session (REST1/REST2), it computes FD and valid scan time after excluding high-motion volumes.
# Subjects with ≥10 minutes of valid data in both sessions are retained.
# Outputs (saved in dir_results):
# - Valid subject lists for LR, RL, and combined encodings (intersection)
# - FD summary per subject/session/encoding
# source(file.path(src_dir,"1_fd_time_filtering.R")) During this step, an FD summary table is generated with the following columns:
subject: HCP subject ID
session: REST1 or REST2
encoding: LR or RL
mean_fd: mean framewise displacement
valid_time_sec: total duration of valid data in seconds
# Read FD summary
fd_summary <- read.csv(file.path(dir_data, "fd_summary.csv"))
# Display the first 4 rows
knitr::kable(head(fd_summary, 4), caption = "First rows of FD summary table")| X | subject | session | encoding | mean_fd | valid_time_sec |
|---|---|---|---|---|---|
| 1 | 100206 | REST1 | LR | 0.1017240 | 858.24 |
| 2 | 100206 | REST2 | LR | 0.1361220 | 858.96 |
| 3 | 100206 | REST1 | RL | 0.0698779 | 864.00 |
| 4 | 100206 | REST2 | RL | 0.0824894 | 863.28 |
As shown above, subject 100206 qualifies for further analysis because each of the four sessions (REST1/REST2 × LR/RL) contains at least 600 seconds of valid data.
The script is currently designed to filter based on valid time only, but it can be easily adapted to apply additional constraints such as maximum mean FD thresholds if desired (e.g., mean_fd < 0.2).
In the final step of subject selection, we balance sex across age groups to reduce potential demographic bias in template estimation.
For the combined list of valid and unrelated subjects,
we:
Subset the HCP unrestricted demographics to include only those subjects.
Split subjects by age group and examine the sex distribution within each group.
If both sexes are present but imbalanced, we randomly remove subjects from the overrepresented group to achieve balance.
Note: If you are not applying the unrelated subject filtering step
(3.2), you can modify the code to subset based on
valid_combined_subjects_FD instead of
valid_combined_subjects_unrelated.
The final list of valid subjects is saved in dir_results
as:
valid_combined_subjects_balanced.csv
valid_combined_subjects_balanced.rds (used in the
template estimation step)
# This script balances sex within each age group for subjects who passed FD and unrelated filtering.
# For each encoding (LR, RL, combined), it samples subjects to equalize the number of males and females per age group, unless an age group includes only one gender (in which case no balancing is applied).
# Uses age and gender information from the unrestricted HCP data.
# Outputs (saved in dir_personal):
# - Sex-balanced subject lists for LR, RL, and combined encodings (as .csv and .rds)
# source(file.path(src_dir,"3_balance_age_sex.R"))In this step, we load and preprocess a group-level cortical
parcellation to be used as the template to estimate the priors in the
next step. Specifically, we use the Yeo 17-network parcellation
(Yeo_17) and perform the following operations:
Simplify the labels by collapsing hemisphere-specific naming and removing subnetwork identifiers, grouping regions by their main network.
Create a new dlabel object that maps each vertex to
its corresponding network.
Mask out the medial wall to exclude it from analysis.
The resulting parcellation is saved as
Yeo17_simplified_mwall.rds.
# This script simplifies the Yeo 17-network parcellation by collapsing region labels and masking out the medial wall.
# It creates a cleaned version of the parcellation suitable for downstream analyses.
# Output:
# - Saved as RDS file in dir_data: "Yeo17_simplified_mwall.rds"
# source(file.path(src_dir,"4_parcellations.R"))[Mandy is here]
We can visualize the Yeo17 networks and their corresponding labels:
# Load libraries
library(ciftiTools)
library(rgl)
rgl::setupKnitr()
# Load the parcellation
yeo17 <- readRDS(file.path(dir_data, "Yeo17_simplified_mwall.rds"))
yeo17 <- add_surf(yeo17)
view_xifti_surface(
xifti = yeo17,
widget = TRUE,
title = "Yeo17 Network Parcellation",
legend_ncol = 6,
legend_fname = "yeo17_legend.png",
)